Prep pGL2-BRI1-GFP seu object with metadata¶

In [1]:
library(tidyverse)
library(Seurat)
library(cowplot)
library(ComplexHeatmap)
library(circlize)
library(GeneOverlap)
library(gprofiler2)
library(ggrepel)
library(ggplot2)
library(muscat)
library(purrr)
library(limma)
library(scran)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.0     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.1     ✔ tibble    3.1.8
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Attaching SeuratObject


Attaching package: ‘cowplot’


The following object is masked from ‘package:lubridate’:

    stamp


Loading required package: grid

========================================
ComplexHeatmap version 2.14.0
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite either one:
- Gu, Z. Complex Heatmap Visualization. iMeta 2022.
- Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
    genomic data. Bioinformatics 2016.


The new InteractiveComplexHeatmap package can directly export static 
complex heatmaps into an interactive Shiny app with zero effort. Have a try!

This message can be suppressed by:
  suppressPackageStartupMessages(library(ComplexHeatmap))
========================================


========================================
circlize version 0.4.15
CRAN page: https://cran.r-project.org/package=circlize
Github page: https://github.com/jokergoo/circlize
Documentation: https://jokergoo.github.io/circlize_book/book/

If you use it in published research, please cite:
Gu, Z. circlize implements and enhances circular visualization
  in R. Bioinformatics 2014.

This message can be suppressed by:
  suppressPackageStartupMessages(library(circlize))
========================================


Loading required package: SingleCellExperiment

Loading required package: SummarizedExperiment

Loading required package: MatrixGenerics

Loading required package: matrixStats


Attaching package: ‘matrixStats’


The following object is masked from ‘package:dplyr’:

    count



Attaching package: ‘MatrixGenerics’


The following objects are masked from ‘package:matrixStats’:

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
    colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
    colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
    colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
    colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
    colWeightedMeans, colWeightedMedians, colWeightedSds,
    colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
    rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
    rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
    rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
    rowWeightedSds, rowWeightedVars


Loading required package: GenomicRanges

Loading required package: stats4

Loading required package: BiocGenerics


Attaching package: ‘BiocGenerics’


The following object is masked from ‘package:limma’:

    plotMA


The following objects are masked from ‘package:lubridate’:

    intersect, setdiff, union


The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union


The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs


The following objects are masked from ‘package:base’:

    anyDuplicated, aperm, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
    get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
    Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
    table, tapply, union, unique, unsplit, which.max, which.min


Loading required package: S4Vectors


Attaching package: ‘S4Vectors’


The following objects are masked from ‘package:lubridate’:

    second, second<-


The following objects are masked from ‘package:dplyr’:

    first, rename


The following object is masked from ‘package:tidyr’:

    expand


The following objects are masked from ‘package:base’:

    expand.grid, I, unname


Loading required package: IRanges


Attaching package: ‘IRanges’


The following object is masked from ‘package:lubridate’:

    %within%


The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice


The following object is masked from ‘package:purrr’:

    reduce


Loading required package: GenomeInfoDb

Loading required package: Biobase

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.



Attaching package: ‘Biobase’


The following object is masked from ‘package:MatrixGenerics’:

    rowMedians


The following objects are masked from ‘package:matrixStats’:

    anyMissing, rowMedians



Attaching package: ‘SummarizedExperiment’


The following object is masked from ‘package:SeuratObject’:

    Assays


The following object is masked from ‘package:Seurat’:

    Assays


Loading required package: scuttle

In [2]:
library(future)
#for 200gb ram 
options(future.globals.maxSize = 200000 * 1024^2)
In [3]:
sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: AlmaLinux 9.3 (Shamrock Pampas Cat)

Matrix products: default
BLAS/LAPACK: /hpc/group/pbenfeylab/tmn23/miniconda3/envs/muscat/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats4    grid      stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] future_1.31.0               scran_1.26.0               
 [3] scuttle_1.8.0               SingleCellExperiment_1.20.0
 [5] SummarizedExperiment_1.28.0 Biobase_2.58.0             
 [7] GenomicRanges_1.50.0        GenomeInfoDb_1.34.8        
 [9] IRanges_2.32.0              S4Vectors_0.36.0           
[11] BiocGenerics_0.44.0         MatrixGenerics_1.10.0      
[13] matrixStats_0.63.0          limma_3.54.0               
[15] muscat_1.12.0               ggrepel_0.9.3              
[17] gprofiler2_0.2.1            GeneOverlap_1.34.0         
[19] circlize_0.4.15             ComplexHeatmap_2.14.0      
[21] cowplot_1.1.1               SeuratObject_4.1.3         
[23] Seurat_4.3.0                lubridate_1.9.2            
[25] forcats_1.0.0               stringr_1.5.0              
[27] dplyr_1.1.0                 purrr_1.0.1                
[29] readr_2.1.4                 tidyr_1.3.0                
[31] tibble_3.1.8                ggplot2_3.4.1              
[33] tidyverse_2.0.0            

loaded via a namespace (and not attached):
  [1] pbdZMQ_0.3-9              scattermore_0.8          
  [3] bit64_4.0.5               irlba_2.3.5.1            
  [5] DelayedArray_0.24.0       data.table_1.14.8        
  [7] KEGGREST_1.38.0           RCurl_1.98-1.10          
  [9] doParallel_1.0.17         generics_0.1.3           
 [11] ScaledMatrix_1.6.0        RhpcBLASctl_0.23-42      
 [13] RSQLite_2.2.20            RANN_2.6.1               
 [15] bit_4.0.5                 tzdb_0.3.0               
 [17] spatstat.data_3.0-0       httpuv_1.6.9             
 [19] viridis_0.6.2             hms_1.1.2                
 [21] evaluate_0.20             promises_1.2.0.1         
 [23] fansi_1.0.4               progress_1.2.2           
 [25] caTools_1.18.2            igraph_1.3.5             
 [27] DBI_1.1.3                 geneplotter_1.76.0       
 [29] htmlwidgets_1.6.1         spatstat.geom_3.0-6      
 [31] ellipsis_0.3.2            backports_1.4.1          
 [33] annotate_1.76.0           aod_1.3.2                
 [35] deldir_1.0-6              sparseMatrixStats_1.10.0 
 [37] vctrs_0.5.2               ROCR_1.0-11              
 [39] abind_1.4-5               cachem_1.0.6             
 [41] withr_2.5.0               progressr_0.13.0         
 [43] sctransform_0.3.5         prettyunits_1.1.1        
 [45] goftest_1.2-3             cluster_2.1.4            
 [47] IRdisplay_1.1             lazyeval_0.2.2           
 [49] crayon_1.5.2              genefilter_1.80.0        
 [51] spatstat.explore_3.0-6    edgeR_3.40.0             
 [53] pkgconfig_2.0.3           nlme_3.1-162             
 [55] vipor_0.4.5               blme_1.0-5               
 [57] rlang_1.0.6               globals_0.16.2           
 [59] lifecycle_1.0.3           miniUI_0.1.1.1           
 [61] rsvd_1.0.5                polyclip_1.10-4          
 [63] lmtest_0.9-40             Matrix_1.5-3             
 [65] IRkernel_1.3.2            boot_1.3-28.1            
 [67] zoo_1.8-11                base64enc_0.1-3          
 [69] beeswarm_0.4.0            ggridges_0.5.4           
 [71] GlobalOptions_0.1.2       png_0.1-8                
 [73] viridisLite_0.4.1         rjson_0.2.21             
 [75] bitops_1.0-7              KernSmooth_2.23-20       
 [77] Biostrings_2.66.0         blob_1.2.3               
 [79] DelayedMatrixStats_1.20.0 shape_1.4.6              
 [81] parallelly_1.34.0         spatstat.random_3.1-3    
 [83] beachmat_2.14.0           scales_1.2.1             
 [85] memoise_2.0.1             magrittr_2.0.3           
 [87] plyr_1.8.8                ica_1.0-3                
 [89] gplots_3.1.3              zlibbioc_1.44.0          
 [91] compiler_4.2.2            dqrng_0.3.0              
 [93] RColorBrewer_1.1-3        clue_0.3-64              
 [95] lme4_1.1-31               DESeq2_1.38.0            
 [97] fitdistrplus_1.1-8        cli_3.6.0                
 [99] XVector_0.38.0            lmerTest_3.1-3           
[101] listenv_0.9.0             patchwork_1.1.2          
[103] pbapply_1.7-0             TMB_1.9.2                
[105] MASS_7.3-58.2             tidyselect_1.2.0         
[107] stringi_1.7.12            BiocSingular_1.14.0      
[109] locfit_1.5-9.7            tools_4.2.2              
[111] timechange_0.2.0          future.apply_1.10.0      
[113] parallel_4.2.2            uuid_1.1-0               
[115] bluster_1.8.0             foreach_1.5.2            
[117] metapod_1.6.0             gridExtra_2.3            
[119] Rtsne_0.16                digest_0.6.31            
[121] shiny_1.7.4               Rcpp_1.0.10              
[123] broom_1.0.3               later_1.3.0              
[125] RcppAnnoy_0.0.20          httr_1.4.4               
[127] AnnotationDbi_1.60.0      Rdpack_2.4               
[129] colorspace_2.1-0          XML_3.99-0.13            
[131] tensor_1.5                reticulate_1.28          
[133] splines_4.2.2             statmod_1.5.0            
[135] uwot_0.1.14               spatstat.utils_3.0-1     
[137] scater_1.26.0             sp_1.6-0                 
[139] plotly_4.10.1             xtable_1.8-4             
[141] jsonlite_1.8.4            nloptr_2.0.3             
[143] R6_2.5.1                  pillar_1.8.1             
[145] htmltools_0.5.4           mime_0.12                
[147] glue_1.6.2                fastmap_1.1.0            
[149] minqa_1.2.5               BiocParallel_1.32.5      
[151] BiocNeighbors_1.16.0      codetools_0.2-19         
[153] utf8_1.2.3                lattice_0.20-45          
[155] spatstat.sparse_3.0-0     numDeriv_2016.8-1.1      
[157] pbkrtest_0.5.2            ggbeeswarm_0.7.1         
[159] leiden_0.4.3              gtools_3.9.4             
[161] survival_3.5-3            glmmTMB_1.1.5            
[163] repr_1.1.6                munsell_0.5.0            
[165] GetoptLong_1.0.5          GenomeInfoDbData_1.2.9   
[167] iterators_1.0.14          variancePartition_1.28.0 
[169] reshape2_1.4.4            gtable_0.3.1             
[171] rbibutils_2.2.13         
In [4]:
rc.integrated <- readRDS("../../CheWei/scRNA-seq/Integrated_Objects/rc.integrated_6S_GL2_lines_seu3_20240423.rds")
In [5]:
rc.integrated
An object of class Seurat 
71741 features across 37615 samples within 3 assays 
Active assay: SCT (24847 features, 0 variable features)
 2 other assays present: RNA, integrated
 3 dimensional reductions calculated: pca, umap, umap_2D
In [6]:
table(rc.integrated$orig.ident)
sc_130 sc_131 sc_133 sc_134 sc_135 sc_137 
  6656   7658   8424   7837   4975   2065 
In [7]:
rc.integrated$geno <- rc.integrated$orig.ident
rc.integrated$geno <- gsub("sc_130","WT",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_134","WT",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_131","bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_135","bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_133","pGL2_BRI1_GFP_bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_137","pGL2_BRI1_GFP_bri1_T",rc.integrated$geno)
In [8]:
rc.integrated$rep <- rc.integrated$orig.ident
rc.integrated$rep <- gsub("sc_130","1",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_134","2",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_131","1",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_135","2",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_133","1",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_137","2",rc.integrated$rep)
In [9]:
table(rc.integrated$orig.ident, rc.integrated$geno)
        
         bri1_T pGL2_BRI1_GFP_bri1_T   WT
  sc_130      0                    0 6656
  sc_131   7658                    0    0
  sc_133      0                 8424    0
  sc_134      0                    0 7837
  sc_135   4975                    0    0
  sc_137      0                 2065    0
In [10]:
feature_names <- read_tsv("../data/features.tsv.gz", col_names = c("AGI", "Name", "Type")) %>%
  select(-Type) %>%
  distinct()
Rows: 32833 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (3): AGI, Name, Type

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
In [11]:
rc.integrated$geno <- factor(rc.integrated$geno, levels=c("WT", "bri1_T", "pGL2_BRI1_GFP_bri1_T"))

Cell and developmental stage metadata¶

  • Developmental stage: time_zone
  • Cell type:cell_type
  • Combination of cell type and developmental stage: time_zone_cell_type
  • Combination of cell type and developmental stage with cell subtypes (not used): time_zone_cell_subtypes
In [12]:
rc.integrated$cell_type <- rc.integrated$celltype.anno.Li.crude
rc.integrated$time_zone <- rc.integrated$time.anno
rc.integrated$time_zone_cell_type <- rc.integrated$time.celltype.anno.Li.crude
In [13]:
table(rc.integrated$orig.ident, rc.integrated$cell_type)
        
         Atrichoblast Columella Cortex Endodermis Lateral Root Cap Pericycle
  sc_130         1468       893    982        520             1078       271
  sc_131         1405       669    800        595             1787       749
  sc_133         1421       605   1452        884             1684       332
  sc_134         1195       588    924        639             2136       665
  sc_135         1040       287    541        525              963       407
  sc_137          351       282    293        236              421       134
        
         Phloem Procambium Trichoblast Xylem
  sc_130     57        190        1021   176
  sc_131    119        331         944   259
  sc_133     54        145        1669   178
  sc_134    123        353         937   277
  sc_135    106        211         707   188
  sc_137     27         68         197    56
In [14]:
table(rc.integrated$time.anno)
         Distal Columella   Distal Lateral Root Cap                Elongation 
                     2943                      2941                     10498 
               Maturation                  Meristem        Proximal Columella 
                     7025                      8801                       359 
Proximal Lateral Root Cap 
                     5048 
In [15]:
# Plot celltype annotation Li
order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem","Protophloem", "Xylem", "Procambium","Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400D3", "#DCD0FF","#5AB953", "#BFEF45", "#008080", "#21B6A8", "#82B6FF", "#0000FF","#E6194B", "#DD77EC", "#9A6324", "#FFE119", "#FF9900", "#FFD4E3", "#9A6324", "#DDAA6F", "#EEEEEE")
rc.integrated$cell_type <- factor(rc.integrated$cell_type, levels = order[sort(match(unique(rc.integrated$cell_type),order))])
color <- palette[sort(match(unique(rc.integrated$cell_type),order))]


options(repr.plot.width=15, repr.plot.height=7)

(BR_cell <- DimPlot(rc.integrated, reduction = "umap", group.by = "cell_type", cols = color, split.by = 'geno', pt.size = 0.75, ncol=3)) 
ggsave("../output/GL2/pGL2_BRI1-GFP_cell_types.pdf", width=13, height=6)
In [16]:
options(repr.plot.width=10, repr.plot.height=6)
FeaturePlot(rc.integrated, features="GL2-BRI1-GFP", split.by = "geno", order=T)
#ggsave("../output/GL2/pGL2_BRI1-GFP_expression.pdf", width=10, height=6)
Warning message in FeaturePlot(rc.integrated, features = "GL2-BRI1-GFP", split.by = "geno", :
“All cells have the same value (0) of GL2-BRI1-GFP.”
In [17]:
options(repr.plot.width=10, repr.plot.height=6)
FeaturePlot(rc.integrated, features="AT1G79840", split.by = "geno", order=T)
#ggsave("../output/GL2/GL2-AT1G79840_expression.pdf", width=10, height=6)
In [18]:
options(repr.plot.width=30, repr.plot.height=7)

(BR_cell <- DimPlot(rc.integrated, reduction = "umap", group.by = "cell_type", cols = color, split.by = 'orig.ident', pt.size = 0.75, ncol=8))
In [19]:
DefaultAssay(rc.integrated) <- "SCT"
FeaturePlot(rc.integrated, features="GL2-BRI1-GFP", split.by = "orig.ident", order=T)
Warning message in FeaturePlot(rc.integrated, features = "GL2-BRI1-GFP", split.by = "orig.ident", :
“All cells have the same value (0) of GL2-BRI1-GFP.”
Warning message in FeaturePlot(rc.integrated, features = "GL2-BRI1-GFP", split.by = "orig.ident", :
“All cells have the same value (0) of GL2-BRI1-GFP.”
In [20]:
saveRDS(rc.integrated, file = "/hpc/group/pbenfeylab/CheWei/scRNA-seq/Integrated_Objects/rc.integrated_6S_GL2_lines_seu4_annotated_20240429.rds")
In [ ]: